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We present improved results for the B and D meson spectrum from lattice QCD including the 
effect of u/d,s and c quarks in the sea. For the B mesons the Highly Improved Staggered Quark 
action is used for the sea and light valence quarks and NonRelativistic QCD for the b quark including 
0(a s ) radiative corrections to many of the Wilson coefficients for the first time. The D mesons 
use the Highly Improved Staggered Quark action for both valence quarks on the same sea. We find 
M Bs -M B = 84(2) MeV, M Bs = 5.366(8) GeV, M Ba = 6.278(9) GeV, M D , = 1.9697(33) GeV, and 
Md s - M D = 101(3) MeV. Our results for the B meson hyperfine splittings are M B * - Mb = 50(3) 
MeV, Mb* — Mb b = 52(3) MeV, in good agreement with existing experimental results. This 
demonstrates that our perturbative improvement of the NRQCD chromo-magnetic coupling works 
for both heavyonium and heavy-light mesons. We predict Mb* —Mb c = 54(3) MeV. We also present 
first results for the radially excited B c states as well as the orbitally excited scalar B* and axial 
vector _B c i mesons. 



I. INTRODUCTION 

Lattice QCD calculations have become an essential 
part of B physics phenomenology [T], providing increas- 
ingly precise determinations of decay constants and mix- 
ing parameters needed, along with experiment, in the 
determination of CKM matrix elements. Since these cal- 
culations can now give stringent constraints on the CKM 
unitarity triangle, currently resulting in tension at a few 
sigma level [2], it is important to check that all system- 
atic errors have been correctly accounted for. With this 
in mind we present a new study of the B-meson spectrum 
that provides a good check of recent improvements that 
have been made in our discretization of the QCD La- 
grangian. The possibility of more B states being found 
at experiments such as LHCb also gives us the oppor- 
tunity for further tests of QCD in the nonpcrturbative 
regime. We emphasise that all parameters for this cal- 
culation, including quark masses and the lattice spacing, 
have already been determined elsewhere [3J making this 
a parameter free test of lattice QCD. 

This test is made possible by the use of NonRelativistic 
QCD (NRQCD) for the b quark, which has the advan- 
tage that the same action can be used for both bottomo- 
nium and B-meson calculations. HPQCD recently com- 
puted the one loop radiative corrections to many of the 
coefficients in the NRQCD action [3J 0] and studied the 
effect of these improvements on the bottomonium spec- 
trum in [3j. Systematic errors were significantly reduced 
in a number of quantities, including the hyperfine split- 
ting, and the first QCD prediction of the D-wave spin 
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splittings was made [5|. This analysis used new gluon 
configurations (ij generated by the MILC collaboration 
with 2+1+1 flavours of HPQCD 's Highly Improved Stag- 
gered Quarks (HISQ) [7] in the sea and including n/a s a 2 
improvements to the gluon action [8]. We use the same 
gluon configurations here. 



For the u, d, s and c valence quarks in our calcula- 
tion we use the same HISQ action as for the sea quarks. 
The advantage of using HISQ is that am q discretisation 
errors are under sufficient control that it can be used 
both for light and for c quarks 7J. Both the NRQCD 
and the HISQ actions are also numerically very cheap 
which means we are able to perform a very high statistics 
calculation to combat the signal to noise ratio problems 
that arise in simulating B-mesons. The same u/d,s,c 
HISQ quark propagators used in the B mesons can also 
be used to calculate the masses of pseudoscalar charmed 
mesons which we also present here. Our results are pre- 
cise enough that it is possible to distinguish the heavy 
quark dependence of splittings such as the Md s — Mo 
and M B - M B - 



We begin by outlining the methods used in our lat- 
tice calculation, which are similar to [3j [9] . The B s , B c 
and B meson masses and the radially excited B c are pre- 
sented in Sec. |III| hyperfine results are given in Sec. [TV 
axial vector and scalar B-mesons are discussed in IV 
Sec. |VI| compares our results to earlier NRQCD-HISQ 
ones on nj =2 + 1 configurations including asqtad sea 
quarks 9 and to calculations using the HISQ action for b 
quarks |10Ullj . Sec. VII gives our conclusions, including 
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an updated spectrum for gold-plated mesons from lattice 
QCD. 
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II. LATTICE CALCULATION 

Our calculation uses five ensembles of gluon configu- 
rations generated by the MILC collaboration [6]. These 
are rif = 2 + 1 + 1 configurations that include the effect 
of light, strange and charm quarks with the HISQ action 
and a Symanzik improved gluon action with coefficients 
correct through 0(a s a 2 , rifa s a 2 ) [8]. The lattice spacing 
values range from a — 0.15fm to a — 0.09fm. The config- 
urations have accurately tuned sea strange quark masses 
and sea light quark masses (m u = m ( i = mi) with ra- 
tios to the strange mass of mi/m s = 0.1 and 0.2, which 
correspond to pions of mass 220-315 MeV. Having sea 
quark masses close to the physical point is particularly 
important for studies of the B meson where chiral ex- 
trapolations make up a substantial portion of the final 
error. 

In Ref. [3] we accurately determined the lattice spac- 
ings using the T(25 — IS") splitting and the decay con- 
stant of the fictitious 7] s particle, a pseudo-scalar ss me- 
son whose valence quarks are not allowed to annihilate 
on the lattice [T2] . Agreement was shown between these 
methods in the continuum limit. In this paper we use the 
T(2S — IS") lattice spacings. The details of each ensem- 
ble, including the sea quark masses and spatial volumes, 
are given in Table [i] All ensembles were fixed to Coulomb 
gauge. 

Light, strange and charm quark propagators were gen- 
erated using the HISQ action, the masses used are given 
in Table [Tl] In Ref. [3J accurate strange quark masses 
were given for each ensemble, tuned from the mass of 
the 7y s meson, which was determined from K and 7r me- 
son masses to be 0.6893(12) GeV. The values of am™ 1 in 
Table [Tl] correspond to these. Mistiming of the strange 
quark mass was a major source of error in Ref. [9] which 
will not be present in this calculation. The light valence 
quarks are taken to have the same masses as in the sea. 

Charm quark masses are tuned by matching the mass 
of the rj c to experiment. The experimental value is shifted 
by 2.6 MeV for missing electromagnetic effects and 2.4 
MeV for not allowing it to annihilate to gluons, giving 
2.985(3) GeV [12]. The e Naik term in the action is not 
negligible for charm quarks and we use the tree level for- 
mula given in |13j . the values appropriate to our masses 
are given in Table [H] 

The velocity of a b quark in a bound state is typi- 
cally very small; v 2 = 0.1 in bottomonium and v 2 varies 
from 0.01 to 0.04 in heavy light systems containing a 
b quark. This makes NRQCD [T3] a suitable effective 
field theory for handling b quarks. It also has a num- 
ber of other advantages. By construction, we are able 
to perform calculations at relatively coarse lattice spac- 
ings since discretisation errors are not set by powers of 
the quark mass as in a relativistic theory. Generation of 
propagators is very fast since in NRQCD they can simply 
be generated by time evolution with a given Hamiltonian. 
The other major benefit is that NRQCD can be used for 
both heavy-heavy and heavy-light mesons. All free pa- 



rameters in this calculation were previously tuned using 
the statistically more precise bottomonium spectrum in 
[3J , meaning that all results here are parameter free tests 
of QCD. 

These advantages come at a price. NRQCD is non- 
renormalisable because operators of dimension greater 
than four are included in the action, rather than being 
evaluated as operator insertions as in HQET. This means 
that the continuum limit a — > cannot be taken. This 
does not mean, however, that physical results cannot be 
extracted. Because NRQCD is an effective theory, con- 
tinuum results can be inferred from fits to calculations in 
its regime of validity, where am^ > 1. We discuss this in 
Sec. HI A 1| As finer lattices become more readily avail- 
able on which am^ < 1, other methods |10j may become 
more appropriate than NRQCD. In the meantime, how- 
ever, NRQCD remains the easiest and best way to access 
the full range of heavy quark physics in lattice QCD. 

The NRQCD Hamiltonian we use is given by [15] : 

aH = aHo + a5H ; 
A( 2 ) 
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Here V is the symmetric lattice derivative and A^ 2 ) and 
A' 4 ) the lattice discretization of the continuum £V D 2 
and J2i Df respectively, am^ is the bare b quark mass. E 
and B are the chromoelectric and chromomagnetic fields 
calculated from an improved clover term |16j . The B and 
E are made anti-hermitian but not explicitly traceless, 
to match the perturbative calculations done using this 
action. 

The coefficients Ci in the action are unity at tree level 
but radiative corrections cause them to depend on ami, 
at higher orders in a s . These were calculated for the 
relevant b quark masses using lattice perturbation theory 
in 3J and the values used in this paper are given in Table 
[TIT A major improvement in this work is the inclusion of 
one loop radiative corrections to C4 [I] which controls the 
hyperfine splitting between the vector and pseudo-scalar 
states. We show in Sec. IIVI that this leads to accurate 
results for 6-light hyperfine splittings in keeping with the 
results of [3] for bottomonium. 

The tuning of the b quark mass on these ensembles 
was discussed in [3J. We use the spin-averaged kinetic 
mass of the T and 77b and take the experimental value to 
which we tune to be 9.445(2) GeV. This allows for elec- 
tromagnetism and 7/b annihilation effects missing from 
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TABLE I. Details of the five gauge ensembles used in this calculation [6]. /3 is the gauge coupling, ox is the lattice spacing 
as determined by the T(2S — IS) splitting in j3], where the three errors are statistics, NRQCD systematics and experiment. 
ami,am s and am c are the sea quark masses, L/a x T/a gives the spatial and temporal extent of the lattices and n c f g is the 
number of configurations in each ensemble. The ensembles 1 and 2 will be referred to in the text as "very coarse", 3 and 4 as 
"coarse" and 5 as "fine". 



Set 


P 


ot (fm) 


ami 


am a 


am c 


L/a X T/a 


n c f g 


1 


5.80 


0.1474(5) (14) (2) 


0.013 


0.065 


0.838 


16x48 


1020 


2 


5.80 


0.1463(3)(14)(2) 


0.0064 


0.064 


0.828 


24x48 


1000 


3 


6.00 


0.1219(2) (9) (2) 


0.0102 


0.0509 


0.635 


24x64 


1052 


4 


6.00 


0.1195(3) (9) (2) 


0.00507 


0.0507 


0.628 


32x64 


1000 


5 


6.30 


0.0884(3) (5) (1) 


0.0074 


0.037 


0.440 


32x96 


1008 



TABLE II. 



The parameters used in the generation of the 



HISQ propagators, nmj are the valence quark masses and 



e N aik is the coefficient of the Naik term in the action. 
4 tNaik is very slightly wrong - it should be -0.224. 



On set 



Set amT 1 




amT 1 


amT 1 




1 0.013 




0.0641 


0.826 


-0.345 


2 0.0064 




0.0636 


0.818 


-0.340 


3 0.01044 




0.0522 


0.645 


-0.235 


4 0.00507 




0.0505 


0.627 


-0.222 


5 0.0074 




0.0364 


0.434 


-0.117 


TABLE III. The coefficients ci, c 5 , 


C4 and C6 


used in the 


NRQCD action [T] c 2 and 


C3 are set to 1.0. 




Set 


Cl 


C5 


c 4 


C6 


very coarse 


1.36 


1.21 


1.22 


1.36 


coarse 


1.31 


1.16 


1.20 


1.31 


fine 


1.21 


1.12 


1.16 


1.21 



our calculation [9] . Note that we no longer have to apply 
a shift for missing charm quarks in the sea The val- 
ues used in this calculation are tuned on that basis and 



given in Table IV along with other parameters. 

The calculation of NRQCD-HISQ two point functions 
with stochastic noise sources uses the method developed 
in [9] to allow spin-information to be added into the 
HISQ propagators so that the correct J PC NRQCD- 
light correlators can be made. Once HISQ propaga- 
tors have been made with a source time-slice of random 
numbers we can no longer apply the 'staggering matrix', 
fl(x) = riu=i it-^ i & t the source to convert them to naive 
quark propagators with spin as would be used in the orig- 
inal method for combining staggered and non-staggered 
quarks |17j . Instead we include the staggering matrix at 
the source of the NRQCD propagators along with the 
same time-slice of random numbers [9] I18j . 

We also use exponentially smeared quark sources, 
which take form exp(— r/a sm ) as a function of radial dis- 
tance, for the NRQCD propagators. These use two dif- 
ferent radial sizes, a sm , on each ensemble as given in Ta- 
ble [IV] Correlators were calculated at 16 time sources on 
each configuration and the calculation was repeated with 



TABLE IV. Parameters used in the NRQCD action, amb 
is the bare b quark mass and uol the Landau link tadpole- 
improvement factor used in the NRQCD action [19] . n c f g gives 
the number of configurations used in each ensemble. 16 time 
sources were used on each configuration. The column a sm 
gives the size parameters of the quark smearing functions, 
which take the form exp(— r/a sm ). a sm kept approximately 
constant in physical units. 



Set 


amb 


UOL 


n t 




1 


3.297 


0.8195 


16 


2.0,4.0 


2 


3.263 


0.82015 


16 




3 


2.66 


0.834 


16 


2.5,5.0 


4 


2.62 


0.8349 


16 




5 


1.91 


0.8525 


16 


3.425,6.85 



the heavy quark propagating in the opposite time direc- 
tion. All correlators on the same ensemble were binned 
to avoid underestimating the errors. Our method also 
requires the calculation of T and r\h correlators to sub- 
tract the unphysical ground state energy of NRQCD, for 
details see [3]. 

B meson energies are extracted from the two- 
point functions using a simultaneous multi-exponential 
Bayesian fit (20] [21] to the form 

JVcxp 

C mcson { h 3, to;t)=J2 Kkbl k e- Ek{t - to) (2) 
fe=i 

- d^dl k ,(-l)^e- E '^-^. 



k' = l 

The priors on the energy splittings E n+ i — E n are 
600(300) MeV and the priors on the ground states are es- 
timated from previous results with a width of 300 MeV. 
The priors on the amplitudes are 0.1(1.0) and the fit in- 
cludes points from some to to it/2, half the temporal 
extent of the lattice, i and j label the different source 
and sink smearing functions used in the correlator. t Q is 
taken from 7 — 8 on the fine ensemble, 6 — 8 on coarse and 
for very coarse the B and B s are fit from to = 4 — 8 but 
the B c fits started at to = 14 in order to obtain an ac- 
ceptable fit. The B, B s and B c are fit separately but all 
vector and pseudo-scalar correlators for each meson are 
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included in the same fit. Scalar and axial vector states 
are obtained from the oscillating terms (i.e. the E' k ,) in 
Eq. ([2]). The oscillating terms correspond to opposite 
parity states made by the time-doubled quark and are 
typically present in meson correlators made from stag- 
gered quarks. 



we do not do this here. As Fig. [T] shows, no significant 
lattice spacing or sea quark mass dependence is visible in 
our results for Mb 3 i att . 

1. Extracting physical results 



III. MESON MASSES 



We begin with results for pseudo-scalar mesons. Hy- 



pcrfinc splittings are discussed in Sec. IV and scalars and 
axial vectors in Sec. IVl 



A. The B s meson 

In NRQCD meson energies have an unphysical energy 
shift and we must consider energy splittings in order to 
compare with experiment. We subtract half the spin av- 
erage aE bb of the T and rj b ground state energies from 
aE Bs 



1,1, 



latt 



From this we can reconstruct Mb s using 



M 



B ,l a tt - A Ss + ^^b.phys 



(3) 



(4) 



9.445(2) is the relevant experimental 



where M bbphys 
value. 

Our results for o.Eb s and aE bb are given in Table [V] 
Our b and s quark masses are well-tuned here. Never- 
theless we allow small adjustments to Ab s to allow for 
mistuning. These are based on previous determinations 
of the linear slope of with appropriate meson mass, 
M bb for b and Mj? for s. In [5], the slope of Ab s against 
M bb was found to be 0.017 using two values of amj, on a 
very coarse ensemble. By comparing our spin averaged 
kinetic masses to the experimental value on each ensem- 
ble, we obtain the shift Aj\/ 6 - that needs to be applied 
to Ab„ to give the value at the correct b quark mass. 
Using two values of am s on set 1, we find that the slope 
of Ab s with is 0.24(4), consistent with previous re- 
sults [HI [TO]- Comparing M rjs on each ensemble to the 
physical value of 0.6893(12) GeV in [3] gives the tuning 
shift Ajw2 . This is significantly smaller in all cases than 
the lattice spacing error in Ab s ■ The error on both shifts 
is taken to be half the shift itself. 

The splittings Ab s before shifts are applied are listed 
in Table |VI| along with the shifts due to mistuning. 
MB s ,iatt is plotted in Fig. [T] The error is dominated 
by that from the lattice spacing uncertainty. This error 
would be reduced if we constructed an energy difference 
which was much smaller, for example subtracting M Vs 
from both sides of Eq. pi). However the resulting quan- 
tity would then be very sensitive to the s quark mass, so 



Extracting continuum results from a lattice NRQCD 
calculation is more complicated than in a relativistic for- 
malism due to the way coefficients scale with the cutoff. 
Usually, one appropriately tunes parameters in the ac- 
tion so that the results are independent of the cutoff up 
to some power of a, and then fits the remaining depen- 
dence. For example in an 0(a) improved action, the 
following form would be used: 

f(a) = /p llys (1 + h(Aa) 2 + k 2 (Aa) 4 + ...) , 

where A sets the scale and logarithmic terms are generally 
ignored as they are not distinguishable from powers. 

Our results here have discretization errors of the above 
form from the light quark and gluon actions. On top of 
this, our NRQCD action will have discretization errors 
that could have a mild unphysical dependence on anib 
over the range of am b values we are using here (1.9-3.3), 
well within the range of validity of NRQCD as an effec- 
tive theory. The am b dependence comes from missing 
radiative corrections to discretisation correction terms, 
those with coefficients C5 and in Eq. ([I]). 0{a s ) cor- 
rections to these coefficients arc included here, so the 
missing terms are 0(a 2 s ) and higher. To allow for this, 
we include dependence of the discretisation errors on amj, 
in our fits, using the form: 
2 

/(«) = /pbys[l + 2j^j(Aa) 2j (l + d jb Sx m + dj bb (Sx m ) 2 )]. 
4=1 

Here we model the am b dependence with a polynomial 
using the parameter Sx m = (am b — 2.7)/1.5 which varies 
from approximately -0.5 to 0.5 across the range of am b we 
use. In this way we obtain physical results just as with 
any other quark formalism and the error budget from 
the fit includes the additional error from the effective 
field theory cutoff dependence. Note that the effect of 
relativistic corrections to the NRQCD action, which are 
physical, cannot be judged from fitting the data and are 
included as a separate error item. 

In practice we find that most quantities in this work 
have very small lattice spacing dependence. The quanti- 
ties which do show some dependence are the B c mass and 
hyperfine splitting where we believe that the discretisa- 
tion errors come mainly from the charm quark. 

The complete fit function for A^, also includes terms 
to allow for sea quark mass dependence. We take a poly- 
nomial in the variables 8x s and 5xi, defined as the differ- 
ence from the correct quark mass TO qjSe a,phys normalised 
by the correct s quark mass 



Sx n 



[,sca mq,sca,phys 
?^s,sca,phys 
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Set 


arrib 


aM b - b 


aE Vb 


aEr 


am s 


aM Vs 


aE Bs 


aA h a vp 


— 


1 + 1 - 


\ 


? 9Q7 


i .iiy^y j 














U . iOZ 1 _L Zr J 


9RQ( 1 7^ 

U. ZiOc7 1 J- 1 y 


2 


3.263 


7.040(8) 


0.21546(3) 


0.26669(5) 


0.0636 


0.51078(8) 


0.61132(26) 


0.03705(47) 


0.285(5) 


0.280(8) 


3 


2.66 


5.761(14) 


0.22040(5) 


0.26394(7) 


0.0522 


0.42351(9) 


0.52385(23) 


0.03177(18) 


0.228(3) 


0.225(5) 


4 


2.62 


5.719(7) 


0.22408(3) 


0.26767(5) 


0.0505 


0.41476(6) 


0.52029(17) 


0.03102(16) 


0.218(6) 


0.222(4) 


5 


1.91 


4.264(11) 


0.21519(2) 


0.24802(2) 


0.0364 


0.30884(11) 


0.41051(17) 


0.02310(14) 


0.164(5) 


0.161(6) 



TABLE V. Results for energies and kinetic masses in lattice units needed for the determination of the mass of the B s meson. 
The second column gives the b quark mass used on each set. The third to fifth columns are the spin average of the T and r\b 
kinetic masses along with the ground state energies, the values for sets 3-5 are taken from [3] and use C4 = 1. It was shown in 
[3] that the spin averaged kinetic mass does not depend strongly on C4 and since aM bb is only used for small tuning adjustments 
this value is sufficient. Column 6 gives the strange quark mass used in each run. Column 7 is the mass of the n s meson at 
the corresponding strange mass, again taken from [3], apart from retuning on sets 1 and 2. The ground state energies of the 
pseudoscalar B s are given in column 8 and the hyperfine splitting A!? yp = E(Bl) — E(B B ) in column 9. Columns 10 and 11 
give the values of mass differences between scalar and pseudoscalar and between axial vector and vector respectively. 



Set 


A Ss (GeV) 


Am, (MeV) 


A M 2 (MeV) 


5.i:i 


8x s 


1 


0.6558(7)(67) 


-1.5 


0.0 


0.17 


0.01 


2 


0.6533(4) (67) 


-0.9 


0.1 


0.06 


0.01 


3 


0.6432(4) (47) 


2.0 


1.2 


0.16 


-0.04 


4 


0.6471(3)(49) 


0.0 


1.4 


0.06 


-0.04 


5 


0.6487(4) (44) 


-1.2 


0.0 


0.16 


0.02 



TABLE VI. Results for Ab s (the mass difference between 
the B 3 meson and the spin average of T and r/i, masses) on 
different ensembles. The two errors are statistics and lattice 
spacing uncertainty. Column 3 and 4 give the shifts in MeV 
that are applied to Ab s to compensate for the mistuning of 
the b and s quarks respectively. Errors are 50% of the value 
given. Columns 5 and 6 give 8x1 and 8x 3 , the fractional mis- 
tuning of the sea quark masses in units of the s quark mass, 
as defined in the text. 



The values of 8x q entering the fits arc given in Table VI 
The values of Sx s are significantly smaller than for the 
Asqtad 2+1 ensembles used before [9] and the 8x1 values 
correspondingly closer to the physical point. 

With this chiral dependence included, the fit function 
becomes: 



dence is typically 1/3 of valence mass dependence 
which would give a slope of 0.2 here. 

• Quadratic sea quark mass dependence terms 
bu, bi s , b ss are smaller by another factor of 0.2, giv- 
ing 0.000(13). 

• The remaining a 4 and am b terms, d2,dj b and dj bb , 
are given a wide prior of 0(1). 

The fit gives A^phys = 0.644(6) GeV and is robust un- 
der changes in the priors and fit function. The 6 MeV 
error can be broken down into contributions from a de- 
pendence, sea quark mass dependence and the error on 
the data points by looking at the variation of the x 2 
[21] . These contributions are listed separately in our fi- 
nal error budget and are dominated by the error on the 
data points, i.e. statistics and lattice spacing uncertainty. 
Since the quark masses are very well tuned, the correc- 
tions for mistuning applied in the previous section pro- 
duce negligible effects. 



2. Systematic errors 



A Bs ( a ^xi,Sx s ) = A Bs! p hys [l 
2 

+ ^ dj(Aa) 23 {l + d jb Sx m + d jbb (8x m ) 2 ) 

+ 2b l 8x l (l + d l {ka) 2 ) 
+ b s 8x s (l + d s (Aa) 2 ) 
+ ±b u {8xi) 2 + 2b ls Sx l 8x s + b ss (8x s ) 2 }. (5) 

We take the prior on A^^phys to be 0.6(2) and we 
take the physical scale to be A = 400 MeV based on the 
typical meson momenta. The other terms and priors are: 

• The quadratic a dependence terms di , di , d s should 
be 0(a s ) or smaller and so have a prior 0.0(3). 

• The leading sea quark mass dependence terms 6; , b s 
have priors 0.00(7) since sea quark mass depen- 



We now describe the remaining sources of systematic 
error that cannot be estimated from the fit. The largest 
of these is the spin independent NRQCD systematic error 
although there is a significant improvement over previous 
work due to the inclusion of radiative corrections. 

Spin independent NRQCD systematics: This error 
can affect both the bottomonium and B s pieces of 
A^phys- For bottomonium, the NRQCD action is cor- 
rect through 0(a s v 4 ) so the largest errors will be a 2 v A 
and v 6 . v 2 effects are of order 500 MeV, so we allow an 
error of 0.3 2 x 0.1 x 500 = 4.5 MeV from missing a s v A 
corrections (compared to 15 MeV in [S].) Similarly, v 6 
terms should be 5 MeV. Adding these in quadrature and 
dividing by two gives 3.4 MeV. For the B s , power count- 
ing is in terms of v = A/m b which is even smaller in a 
heavy-light meson and missing spin independent correc- 
tions are negligible. 
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Spin dependent NRQCD systematica: Since the bot- 
tomonium energies are spin averaged, the only contribu- 
tion from spin dependent terms is to the B s mass. With 
the one loop corrections to C4, the dominant error comes 
from radiative corrections to the o.B term and missing 
(A/mb) 2 terms. We take the error to be 30^/4 times the 
hyperfine splitting B* — B s which gives 3 MeV. 

Electromagnetism: The effects of missing electromag- 
netism were estimated in [9] and give a 0.1 MeV error in 
the B s . 

Finite volume effects: Chiral perturbation theory and 
studies of the wave functions of heavy mesons show that 
finite volume errors are negligible for the ensembles used 
here. 

M r}s and M bb : The uncertainty in the n s mass and 
the error from the electromagnetic and annihilation cor- 
rections in M bb also feed into the total error. M„ s has 
an error of 1.2 MeV which using the slope of 0.24 vs M 2 
gives an error of 0.4 MeV to be added to A Bs . The error 
in the adjusted value of M bb = 9.445(2) GeV has negli- 
gible effect on A Bs , but when reconstructing M Bs this 
leads to a 1 MeV error. The error in M bb comes entirely 
from electromagnetism/annihilation as the experimental 
error is negligible. 



Error 


M Bs 


M Bc ,hh 


M Bc ,hs 


Stats/tuning/uncty in a 


4.8 


2.1 


9.5 


Lattice spacing dependence 


2.2 


0.6 


3.5 


mt dependence 


2.8 


2.9 


3.5 


m qtSea dependence 


1.4 


1.4 


5.0 


spin-ind. NRQCD systs. 


3.4 


5.3 


2.3 


spin-dep. NRQCD systs. 


3.0 


3.0 


0.0 


uncty in M Vs 


0.3 




0.7 


em, annihiln in bb 


1.0 


1.0* 


0.0 


em, annihiln in cc 




1.5* 


0.2 


em effects in B 3 or B c 


0.1 


1.0* 


1.0 


em effects in D s 






1.0 


finite volume 


0.0 


0.0 


0.0 


Total (MeV) 


7.7 


8.0 


12.2 



TABLE VII. Full error budget for B s and B c meson masses in 
MeV. The source of each error is described in the text and the 
total error is obtained by adding in quadrature. Starred errors 
are correlated and are added linearly before being squared. 

The systematic errors are summarised in the error bud- 
get in Table |VII| When added in quadrature the total 
systematic error is 4.7 MeV giving a final value of 

M Bs = 5.366(6) (5) GeV, 

which should be compared with the current PDG value 
of 5.3668(2) MeV [22]. This is the best result for this 
quantity from lattice QCD so far. There is a noticable 
improvement over the systematic errors in Ref. [9] but the 
lattice spacing uncertainty remains similar. The results 
are plotted in Fig. [I] 
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FIG. 1. Results for the B 3 meson mass for each ensemble 
against the lattice spacing. The darker shaded band shows 
the 6 MeV error from the fit and the light band includes 
the estimate of systematic errors. Error bars on the data 
points are uncorrelated and come from statistics, quark mass 
mistuning and uncertainty in the lattice spacing. The data 
points have been corrected for mistuning of valence quark 
masses and missing electromagnetism/^j, annihilation effects. 
The experimental value is included in black for comparison. 



B. The D s meson 

Our method for calculating the mass of the D s meson 
closely follows that of [13] . The previous study on MILC 
2+1 AsqTad ensembles included 5 values of the lattice 
spacing down to 0.045 fm and found Md s = 1.9691(32) 
GeV. Here we have only 3 lattice spacings at the coarser 
end of the range so our result will suffer from a larger 
error from the continuum extrapolation. Some other 
systematic errors are smaller, however and our results 
provide an interesting comparison with those in the B 
spectrum. 

To determine Md s , we calculate the splitting Mn, — 
M Vc /2 which has several advantages over determining the 
mass directly. Since the splitting is much smaller than 
the mass, the same relative scale uncertainty translates 
into a much smaller absolute error on the splitting. It was 
shown in |13) that the c quark mass dependence of the 
splitting is small which leads to reduced tuning errors, 
particularly on the coarsest ensembles where discretisa- 
tion errors are large. Finally, the splitting allows for a 
direct comparison with Mb s — M Vb /2 which must be used 
in the NRQCD case due to the unphysical energy shift. 
The r\ c is used rather than the spin averaged cc state 
simply because a staggered vector meson would require 
additional propagators to be generated. M Bg — M Vb /2 
has a slightly increased systematic error over our pre- 
ferred A Bg (eq. Q). 

Like the NRQCD B s correlators, the D s fit function 
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Set 


am c 


am a 


ami 


o,E Vc 


clEd s 


aE D 


aEo B — aEo 


aEo e — aE Vc /2 






1 




DfS41 


o ni ^ 


9 99 c iOR('7 s l 


1 4S79Q('^n 1 l 

± .*±0 i £iiy you J 






c !747 c i('9Q N | 




n ^ 

U . o 


2 


0.818 


0.0636 


0.0064 


2.21032(4) 


1.47559(20) 


1.41258(68) 


0.06300(69) 


0.37043(19) 


0.0 


0.2 


3 


0.645 


0.0522 


0.01044 


1.83967(5) 


1.21934(14) 


1.17112(53) 


0.04822(47) 


0.29950(13) 


0.7 


0.4 


4 


0.627 


0.0505 


0.00507 


1.80351(3) 


1.19554(8) 


1.14112(61) 


0.05442(58) 


0.29379(8) 


0.8 


0.4 


5 


0.434 


0.0364 


0.0074 


1.33307(4) 


0.88212(9) 


0.84682(26) 


0.03530(21) 


0.21559(10) 


0.0 


0.4 



TABLE VIII. Results from charmed meson fits in lattice units. Columns 2-4 give the HISQ quark masses used in the run, 
columns 5-7 give the rj c , D s and D energies with statistical errors only. Columns 8 and 9 give splittings that we use in our fits. 
They have reduced errors over the naive subtraction of earlier columns because correlations are taken into account. Columns 
10 and 11 give the shifts in MeV that are applied to Ad s due to mistuning of the s and c quarks. 



includes oscillating terms coming from the states related 
by parity and, being relativistic, also includes cosh time- 
dependence: 



0.7 



e -E kt + e -E k (T-t) 



Cmcson(i) — ^ ] a k (t 
fc=l 

-i-iy £ ^( e -^ t + e -^( T -*))(6) 

k' = l 

As for the B s fits, the priors on the energy splittings 
E n+ i — E n are taken to be approximately 600(300) MeV 
and the prior on the ground state is 1.9 GeV with a 300 
MeV width. Similarly the prior splitting between the 
ground state and first oscillating state is 600(300) MeV. 
Fits with iV exp = 5 are typically used as the results are 
stable by this point. D S ,D and r\ c correlators are fit 
simultaneously on each ensemble to include the correla- 
tions in the splittings. 

Before performing a continuum extrapolation we must 
correct for mistuning of the valence quark masses. The s 
and c quark mass dependence of Mjj b — M Va /2 was stud- 
ied in detail in Ref. |13j by fitting the splitting as a func- 
tion of M 2 U and M Vc . The dependence is linear over the 
range of values used with a slope of 0.20(1) against M 2 
and 0.05 against M Vc . Although this data used AsqTad 
sea quarks the corrections are small and since all shifts 
are applied with a 50% error, any difference between the 
slope for HISQ sea quarks will be negligi ble. T he shifts 
applied to Mo, ~M Va /2 are listed in Table VIII Another 
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FIG. 2. Plot of the splittings M Ds - M Vc /2 and M Bs - M Vb /2 
against the square of the lattice spacing. The dark grey band 
is the statistical error and the lighter band gives the full com- 
bined statistical and systematic error. The errorbars include 
statistics, scale and tuning only, correlated errors are not 
shown. The lattice results are adjusted for missing electro- 
magnetic and annihilation effects. 



to have prior 0.5(2) GeV. The discretisation terms have 
priors 0.0(2) except for di which is 0.00(6) since tree- 



advantage of using M Ds - M V J2 is that thccrror from level a 2 ^ errors have been removed in the HISQ action 



the lattice spacing is a third of the naive value. Changing 
the lattice spacing requires m s and m c to be retuned, the 
effect of which partially cancels in the splitting. 

The results at different lattice spacings and light quark 
masses are fit to the same function as in Ref. [13] 

4 

A Ds (a,6xi,5x s ) = A_D S! p hys [l + ^d j (m c a) 2: ' 

+ 2bi5xi(l + di(m c a) 2 ) + b s Sx s (l + d s (m c a) ) 

+ Ab u {5xi) 2 + 2bi s 6xi6x s + b ss (Sx s ) 2 }. (7) 

The same prior values as for the B s are used for the sea 
quark mass dependence and the splitting itself is taken 



Discretisation errors are set by the scale A = m c since 
the dominant error will come from the charm quarks. 

The result of the fit, 0.4808(28) GeV, is plotted in 
Fig. [2] along with the retuned data on each ensemble. 
Also included for comparison is the corresponding split- 
ting in the B meson spectrum Mg a — M Vb /2 from our 
results. There is a significant difference in the two split- 
tings, largely driven by the stronger binding of heavy- 
onium as the heavy quark mass is increased. The ex- 
perimental difference between c and b is well reproduced 
by our results here. The complete dependence on heavy 
quark mass is mapped out in [lOl [11] . The lighter shaded 
band in Fig. [2] includes the systematic errors which are 
discussed in the next section. 



Systematic errors 



eraged bottomonium ground state energy. 



The error arising from statistical/scale, lattice spacing 
dependence and sea quark mass effects is estimated from 
the fit as above, the remaining systematic errors that 
cannot be found in this way are the following: 

Electromagnetism: Electromagnetic effects in the D s 
were estimated in Ref. [13] where the shift was 1-3(7) 
MeV, assuming a 50% error. 

M Vs : The uncertainty in the mass of the rj s meson, 
used for tuning to the correct s quark mass, feeds into the 
error. Using the slope of 0.2, the mass M Vs — 0.6893(12) 
GeV results in an error of 0.2 x 2 x 1.2 x 0.69 - 0.3 MeV. 

M rja : When the D s mass is reconstructed from the 
splitting, we must include the error from M Vc = 2.985(3) 
that comes from our estimate of electromagnetic and an- 
nihilation effects in the r\ c as well as experimental errors. 
This gives a 1.5 MeV error in Mr) s 

Lattice spacing systematics: Systematic errors in the 
determination of the lattice spacing are included in the 
scale error. 



The error budget for Mo s is given in Table IX and our 
final result is 

M Ds = 1.9697(28)(17) GeV 

where the two errors are fitting/scale/tuning and system- 
atics and lead to a combined error of 3.3 MeV. In fact our 
final error is not significantly worse than in [13] because 
an increased lattice spacing extrapolation error is offset 
by the accurate physical value for M rjs . The current ex- 
perimental result for M Ds is 1. 9685(3)" GeV 



Error 


M Da 


Stats/tuning/uncty in a 


2.2 


Lattice spacing dependence 


1.6 


m q , sea dependence 


0.7 


uncty in M Vs 


0.7 


em effects in D a 


0.7 


em, annihiln effects in n c 


1.5 


finite volume 


0.0 


Total (MeV) 


3.3 



TABLE IX. Full error budget for M Ds in MeV. The different 
errors are described in more detail in the text. 



C. The B c meson 

In Ref. [5] , two different methods of reconstructing the 
B c mass were used: the "heavy-heavy" (or hh) subtrac- 
tion method and the "heavy-strange" (or hs) subtraction. 
In the hh method, half the mass of the r\ c is subtracted 
from the lattice value of E Bc in addition to the spin av- 



Mb c = (aE Bc - ^(aE b - b + aM Vc ] 

+ \ ( M bb^ys + Mfe.pby) • 



latt 



(8) 



This has two advantages, firstly it makes the splitting a 
very small value which results in a reduced error from 
the uncertainty in the lattice spacing, and secondly it 
reduces mistuning errors since to a good approximation 
E Ba and M r , c depend linearly on the charm quark mass. 
The second method, hs, uses the B s and D s energies to 
remove the unphysical energy shift from NRQCD: 



M Bc = {aE Bc - (aE Bs + aM D J) 

latt ' 

,phys + -AfD 3 ,phys) • 



(9) 



The D s and r\ c masses are calculated using HISQ for both 
the c and s valence quarks with the parameters given in 
Table [TTJ. The hh and hs methods have different system- 
atic errors and give two independent results to check con- 
sistency. Previously [5] the hh and hs methods resulted 
in total errors in the B c mass of 10 MeV and 19 MeV 
respectively, using NRQCD b quarks. 

Table|X]gives the energies of the B c , D s and r\ c required 
for the two methods. The B s energies are those given in 
Table IVl 



1. Heavy-heavy method 

We begin with the hh method, values for A Bct hh are 
listed in Table |XI[ As for the B s we need to correct 
A-B a ,hh for small mistunings in the quark masses. In [5] 
the slope with respect to M bb was 0.014 (agreeing with 
that from using HISQ b quarks in |TT]) which gives us the 
shifts AM blt hh given in Table 



XI 



The shifts are around 1 
MeV which is comparable to, or slightly larger than, the 
lattice spacing uncertainty. Since the slope is a physical 
dependence rather than a lattice artefact, for the charm 
quark we use the slope against M rja of -0.035 found in 
|11) . This was based on more data and on finer lattice 
spacings than the smaller value in [9]. From this we ob- 
tain the shifts, Am j,/, , in Table IXIl These shifts are 



. Am, c ,Mi in Table _ 
negligible compared to the lattice spacing errors. Again, 
the errors on the shifts are taken to be 50% of the shift. 
As for the B s , once retuning is taken into account the 
actual scale error on the splitting is less than the naive 
value, in this case ranging from 0.5-0.7 of the naive value. 
We take 0.7 times the a error on all ensembles. 

The data are fit to a similar form to that of A Bs but 
with a few changes. Since A Bc j l / l has such a small value, 
the scale, cutoff and sea quark mass dependence are in- 
cluded additively rather than multiplicatively to allow 
them a larger range. We give them instead an overall 
coefficient of 0.4 GeV. We also expect the discretisation 
errors to be dominated by the charm quark so m c ~ 1 
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Set arrib am c am s aE Vc aEo s o,Eb c aA^ yp a A° B ° a A B 1 aE B i — aE Bc aE B ,i — aE B * 

1 3.297 0.826 0.0641 2.22508(7) 1.48729(30) 1.30409(14) 0.03659(17) 0.256(87) 0.212(73) - 

2 3.263 0.818 0.0636 2.21032(4) 1.47559(20) 1.29702(10) 0.03658(13) 0.335(27) 0.327(39) - 

3 2.66 0.645 0.0522 1.83967(5) 1.21934(14) 1.08866(5) 0.03140(3) 0.241(17) 0.250(9) 0.618(27) 0.605(19) 

4 2.62 0.627 0.0505 1.80351(3) 1.19554(8) 1.07252(4) 0.03137(2) 0.252(5) 0.242(6) 0.637(15) 0.625(13) 

5 1.91 0.434 0.0364 1.33307(4) 0.88212(9) 0.81480(3) 0.02470(2) 0.190(2) 0.184(2) 0.616(8) 0.591(7) 



TABLE X. Parameters and results of the B c meson mass calculations. The first three columns give the bottom, charm and 
strange quark masses used in the runs. aM Va , clMd b are the masses of the pseudoscalar mesons generated with the same HISQ 
propagators used for the B c and B s . aE Ba and aA^ p are the ground state energy and hyperfine splitting on each ensemble. 
Columns 9 and 10 give the splittings with the parity partner states discussed in Sec. [V] The final two columns give radial 
excitation energies for the B c and B*. 



Set 


A Bc , hh (GeV) 


&M hh ,hh (MeV) 


A Mflr ,hh (MeV) 


1 


0.0882(2)(9) 


-1.2 


-0.2 


2 


0.0876(2)(9) 


-0.7 


-0.1 


3 


0.0685(1)(5) 


1.7 


-0.2 


4 


0.06999(7) (52) 


0.02 


-0.2 


5 


0.06331(9)(43) 


-1.0 


-0.3 



TABLE XL Table shows the energy splittings Aj^ti where 
the two errors are statistics and lattice spacing uncertainty. 
The second and third columns are the shifts in MeV applied 
to AB c ,hh to adjust for b and c quark mass mistuning respec- 
tively. 



GeV is used instead of A to set their scale. Our fit form 
is then: 

A Bc ,fefe(a, Sx h Sx s ) = A Bc Mh, P hy S + (10) 

4 

0.4[^ dj{m c a) 2] { 1 +d. jb 8x m + d jbb (Sx m ) 2 ) 

3 = 1 

+2biSxi( 1 +di(m c a) 2 + d H (m c a) 4 ) 
+2b s Sx s ( 1 +d B (m c a) 2 + d ss {m c a) A ) 
+ib H (Sx l ) 2 + 2b ls 5x l 5x s +b ss (5x s ) 2 }. 

We take the prior on A,B ci h/i, P hys to be 0.05(5). The pri- 
ors for the fit terms are the same as for the B s case with 
the additional dj,dj b ,dj bb terms having priors of 0(1). 
The fit gives A Bo ^ hhiPhys = 0.06131(39) (fit error only), 
the systematic errors that must be included when recon- 
structing Mb c are the following: 

Spin independent NRQCD systematics: The effect of 
missing terms in the action on M bb is the same as dis- 
cussed previously, but since the b quark velocity in the 
B c is half that in bottomonium we expect partial cancel- 
lation of the a 2 v A errors in As c: hh- We take 1.7 MeV, 
which is half the value for the B s case. The v 6 terms are 
not expected to cancel and results in the same 5 MeV 
giving a total of 2.6 MeV when added in quadrature and 
halved. 

Spin dependent NRQCD systematics: As for the B s , 
we take the error to be 3a^ / 4 times the hyperfine splitting 
in the B c system giving 3 MeV. 
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FIG. 3. Results for the B c meson mass for each ensemble plot- 
ted against the lattice spacing. The errors on the data points 
include statistics, fitting and lattice spacing uncertainty and 
are adjusted for electromagnetic corrections and mistuning of 
quark masses. The dark shaded band gives the error com- 
ing from the fit and the lighter band includes all systematic 
errors discussed in the text. The black circle is the current 
experimental value. 



Electromagnetism: Electromagnetic effects are not 
negligible in the B c and the required shift was estimated 
in [S] to be +2(1) MeV. 

Finite volume effects: Chiral perturbation theory and 
studies of the wave functions of heavy mesons show that 
finite volume errors are negligible for the ensembles used 
here. 

M Vc and M bb : Since the slopes of A.B a ,hh against 
these meson masses are very small, the uncertainty in 



M nc and M bb does not require an additional error to be 
included. However the errors will appear when Mg c is 
reconstructed. These errors come from corrections due to 
electromagnetism and annihilation effects and are corre- 
lated since the same method was used to estimate these 
shifts. Taking half the error on these shifts gives 1 MeV 
for M bb and 1.5 MeV for M Vc . These are added linearly 
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along with the 1 MeV for electromagnetic effects in the 
B c described above before being added in quadrature to 
the other errors. The correlated errors are marked with 
a * in Table ELD 

Taking all of these systematic errors into account, our 
value for the B c mass using the hh method is 

M Bc = 6.278(4) (8) GcV. 

The fit result is plotted in Fig. [3] along with the retuned 
data points for each ensemble. Table |VII| gives the con- 
tribution to the final error of statistics, tuning, scale un- 
certainty and quark mass dependence. Adding the sta- 
tistical and systematic errors in quadrature gives a total 
error of 9 MeV which is shown as the lighter shaded band 
in Fig. 3j The current experimental value is 6.277(6) 
GeV m . 




0.025 



2. Heavy-strange method 

The hs method requires tuning adjustments for the 6, c 
and s quark masses. Ref. [§] found strong dependence 
on the s quark but very small dependence on the b and c 
masses. The slope against M 2 s is 0.41, the slope against 
M b i is 0.005 and against M„ c is 0.07. These slopes agree 
with the results in [11] . The resulting shifts are given in 
Table XII along with the energy splittings A-B e ,hs- The 



biggest shifts are those for mistuning of the s quark on 
the coarse lattices, but even there the shifts are smaller 
than the lattice spacing uncertainty. 



Set 


& Bc ,hs (GeV) 


&M bi ,hs 


Ajm, c ,Ju 

(MeV) 


Am, s ,hs 


1 

2 


-1.069(1)(10) 
-1.065(1)(10) 


-0.4 
-0.25 


0.4 
0.3 


0.0 
0.2 


3 
4 


-1.059(1)(7) 
-1.062(1)(8) 


0.6 
0.0 


0.5 
0.5 


2.1 
2.5 


5 


-1.067(1)(7) 


-0.4 


0.7 


-0.1 



TABLE XII. Results for the hs splitting A Bc ,hs in GeV where 
the two errors are statistical and scale uncertainty. Columns 
3-5 are the shifts in MeV applied for mistuning of the b, c and 
s quark respectively in MeV. 

The fit function is the same as the hh case but with 
the dependences included multiplicativcly 

A Bc ,f ls (a, 5x u 8x s ) = A Bc ,hs, P hy S [l + (11) 

4 

^d,-(m c a) 2j ( 1 +d jb 8x m + d jbb (Sx m ) 2 ) 

3 = 1 

+2biSxi( 1 +di(m c a) 2 + du(m c a) 4 ) 
+b s Sx s ( 1 +d s (m c a) 2 + d ss (m c a) 4 ) 
+Ab u (8x l ) 2 + 2b ls Sx l 8x s + b ss (Sx s ) 2 ] . 

The prior on ^B ,hs,phys is -1-0(2), all other priors are the 
same as for the hh method. The fit result is A.B a ,hs,phys = 



FIG. 4. Results for the B c meson mass for each ensemble 
plotted against the lattice spacing using the hs method. The 
data points are adjusted for missing electromagnetism and 
mistuning of quark masses. The dark grey band is the statis- 
tical error on the fit result and the systematic error is shown 
in light grey and is barely visible on this scale. 



— 1.071(12) where the error is from the fit only. The 
systematic errors are listed below: 

Spin independent NRQCD systematics: There will be 
no cancellation as in the hh case so spin independent 
systematic errors in the B c will be of order a 2 v 4 . Based 
on the b quark velocities in each meson, this should be 
half as big for the B c as for M bb estimated earlier, giving 
2.3 MeV. These missing terms also enter the B s mass but 
are negligible, along with v e terms in both mesons. 

Spin dependent NRQCD systematics: The a ■ B term 
in the action will affect the B s and B c in a similar way so 
errors from unknown a 2 terms in C4 should be negligible. 

Electromagnetism: As in the hh method, there is a 
shift of +2(1) MeV for the B c but we must also include 
a shift of -1.3(7) due to the subtraction of the D s mass. 

M Vs , M Vc and M bb : The errors in the retuning com- 
ing from the M nc and M bb are negligible due to the small 
slopes, but the 1.2 MeV error in the M Vs results in a 0.7 
MeV error in Ab c ,hs- We also need to include the error 
in the reference B s and D s masses which is dominated 
by our estimates of electromagnetic corrections. The D s 
has an error of 0.7 MeV and it is negligible for the B s . 
Including the errors in a correlated way is not necessary 
here as only the electromagnetic shift in the B c is not 
negligible. 

Our final answer for the B c mass with the hs method 



M Bc = 6.264(12) (3) GeV 

where the error is dominated by statistics. This is in good 
agreement, but not quite as accurate, as our result from 
the hh method. This mass is shown in Fig. [4] along with 
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FIG. 5. Energies of the 0~ state in the B c fit against the 
number of exponentials in the fit on the fine ensemble, set 
5. This shows how the value and the error on the first radial 
excitation energy, ai?2, stabilize as the number of exponentials 
is increased. 



the retuned data points on each ensemble, both corrected 
for missing electromagnetic effects described above. 



3. Radially excited states 

Our B c meson correlators fits are accurate enough, and 
include multiple smearings to improve projection on the 
ground state, that there is a good signal for the first 
radially excited states, the B' c and B* . Fig. [5] shows 
how our fit results for these states converge. Unlike the 
B and the B s these states are well below the threshold for 
strong decay, in this case into B, D, so that the states can 
be extracted unambiguously from a lattice calculation 
involving only operators that overlap onto single hadron 
states. The splittings from the ground state are listed in 
Table [X] We only have a signal for the coarse and fine 
ensembles since the starting time in the very coarse fits 
was set too high to extract excited states reliably. 

The splittings between the first radial excitation and 
the ground state are fit to the same form as the hs method 
with a prior on the physical value of 0.5(5). 



in Eq. 11 



Radial splittings are typically very insensitive to quark 
masses so we do not apply any shifts for mistuning. Such 
a shift would be dwarfed by the large statistical errors 
on the splittings. The only significant systematic er- 
ror comes from missing radiative corrections to the spin- 
dependent terms in the action. Since the ground state 
and radially excited state will be affected by this error 
in a similar way we take half the error applied in Sec. 
|IIIC1| giving 1.2 MeV. 

The results from the fits are shown in Figs. [6] and [7J 
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FIG. 6. Results for the splitting M B i — Mb c on the fine and 
coarse ensembles along with the result of the fit. 



our results are: 



Mr 



M Bc = 616(19) sta . t (l) ay8 , MeV 

M B , = 591(18) atat (l) Byst MeV, (12) 



where the error comes almost entirely from statis- 
tics/fitting. The size of these splittings means that we 
expect the mesons to be sufficiently below threshold for 
strong decay into a BD pair to be treated as gold-plated. 
We are unable to resolve the excited hyperfine splitting. 

The radial excitation energies for the B c can be com- 
pared to those for r\ c and r\h- For the rjb recent Belle 
results [23] give 0.597 GeV and for the r\ c the experimen- 
tal average is 0.658 GeV [55]. Our B c result is between 
these two, as might be expected. For the T the experi- 
mental 25 — IS splitting is 0.563 GeV and for the J/ip, 
0.589 GeV 22 . Our B* result agrees reasonably with 
either of these. 



D. The B meson 

We extract the mass of the B meson using the splitting 
Ab = Mb s — Mb in which NRQCD systematics should 
cancel. The mass of the B can then be reconstructed 
using our determination of Mb s in Sec. Ill A Results 
for the lattice energy splittings aEB B — clEb are given in 



Table [Xnl] along with the values of M w on each ensemble 
needed for extrapolation in the light quark mass. The 
large correlation matrix meant that the correlators for 
each meson had to be fit separately but the statistical 
errors are a significant improvement over those in [5]. 

Heavy meson chiral perturbation theory (HM^PT) is 
used for the chiral fits. We use the 1-loop formulas given 
by Jenkins in [24] including heavy quark spin symmetry 
breaking terms at order l/m<g, and up to 0(M 3 ) in the 
light mesons masses. Using the same notation as |24j . 
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Set am t 


am 3 


ami 


aAf, 


aE(B s ) - aE(B) 


aA^ p 




aA^" 1 " 


1 S 9Q7 

2 3.263 

3 2.66 

4 2.62 

5 1.91 


fl 0R41 

U.UU41 

0.0636 
0.0522 
0.0505 
0.0364 


0.0064 
0.01044 
0.00507 
0.0074 


0.16615(7) 
0.19153(9) 
0.13413(5) 
0.14070(9) 


0.05821(110) 
0.04288(63) 
0.04705(60) 
0.03134(78) 


0.0377(9) 
0.0324(4) 
0.0309(4) 
0.0212(11) 


0.207(25) 
0.193(13) 
0.200(4) 
0.159(8) 


9^1 (90') 

0.150(57) 
0.192(15) 
0.207(4) 
0.158(7) 



TABLE XIII. Results in lattice units needed to determine the B meson mass. The first three columns give the b, s and I 
valence quark masses used in the runs. aM n is the pion mass calculated in [3] to be used in the chiral fits. aE(B s ) — aE(Bi) 
is the splitting between the B 8 and B and fflA^ yp is the B hyperfine splitting. The final two columns give the splittings with 
the parity partner states discussed in Sec. [V] 
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FIG. 7. Results for the splitting Mb*' - Mb* on the fine 
and coarse ensembles with the result of the fit. 



the full 5(7(3) formula is 
M Bs - M Bd = -|(2a + 2A^)(m s - m,) 



(13) 



2 

g 7r 

"AT 
3g 2 A 
4A 2 



-Mi - 2M 



j3 



K 



-Mi 



-h{Ml) + l{Ml) + \l{M*) 



where A x = 4^/^ is the chiral scale, a and A^) are 
coefficients of the tree level terms, g is the BB*tt cou- 
pling and A is the coefficient of the term in the effective 
Lagrangian that gives rise to the heavy meson hyperfine 
splitting. The chiral logarithms are given by 



/( If 2 ) l/~ (hx^- + 6 FV (ML) 



including the finite volume correction 
S FV (M L) 



4 ^ Ki{\n\ML) 
M~L ^ \H\ : 



(14) 



(15) 



where K\ is a modified Bessel function and the sum is 
over spatial vectors with components rij £ Z. The finite 



volume corrections shift the pion chiral logarithms by a 
few percent on some ensembles but have a completely 
negligible effect on the fit result. We use the kaon and 
pion masses calculated in [3], and use the tree level rela- 
tion to change M 2 into a combination of M\ and M 2 . 
Quark masses are converted to meson masses using tree 
level relations. 

Our central result uses the reduced SU(2) version of 
the formula: 



C- 



3 (2a- 
4 V 



3g 2 A 
4A 2 



4-2AW) 



mi 



9 2 ir 
A 2 



(16) 



for some constant C. We also perform the fits using the 
5(7(3) formula as a check of systematic errors. Since we 
have a single pion mass for each ensemble and the sea 
strange quark masses are well tuned, partial quenching 
will be a small effect and we use only the full QCD form. 
Staggered quark and other discretisation effects could be 
more significant, however, so the fit function is multiplied 
by 



(1.0 + di(Aa) 2 + d 2 (Aa) 4 ) 



(17) 



at a scale of A = 0.4 GeV. 

We take the prior on g to be 0.5(5) which based on 
several recent lattice calculations [26H29] with a wide er- 
ror covering all of the central values. Our results are 
not sufficient to constrain g, so we test the dependence 
of the final answer on this prior by varying its width. 
While this affects the shape of the curve, the result at 
the physical point does not change significantly since we 
have sufficiently light pion masses. The prior on the tree 
level quark mass term is taken to be 0.5(5) and the pri- 
ors on the discretisation terms di,d,2 are 0.0(5) and 0(1) 
respectively. 

The result of the SU(2) fit is M Bs - M Bd = 85(2) 
MeV when evaluated at a = and at the physical mass 
of the 7T° meson of 0.135 GeV. The fit is shown in Fig. [8] 
and gives a result around la below experiment. To check 
the reliability of the fit, the results of several different fit 
functions are plotted in Fig. [9j This includes the 1-loop 
5(7(2) case, 5(7(2) with different prior widths on g, the 
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FIG. 8. Chiral fit for M Bs - M Bl against the ratio M%/M% s 
for the SU(2) fit function with discretisation effects. The 
grey band shows the chiral fit evaluated at a = with no 
other systematic errors. The plot includes the shift due to 
electromagnetic effects missing in lattice QCD. 
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FIG. 9. Comparison of different chiral fit functions: SU(2) 
HMxPT with and without discretisation corrections; leading 
order analytic terms with and without discretisation terms; 
SU(2) HMxPT with a tighter prior on g of 0.5(1); SU(3) 
HM^PT with and without discretisation terms. Only the 
error from the chiral fit is shown. 



SU(3) case and just the tree level terms with discreti- 
sation effects added in each case. Good x 2 values and 
consistent results are obtained for all fits. 

We now need to consider the effect of electromagnetism 
on As phys - Since our light quark masses are degenerate 
we do not distinguish between the Bd and B u mesons 
but compare to the average M Bl = (M B ± + Mb q )/2. 
Electromagnetism will affect the two states differently 
since the B u is charged. In [3] the shift was estimated to 
be +2(1) MeV for the B u whereas the shift was negligible 



FIG. 10. Plot of the chiral extrapolation for Md b — Mb with 
Mb s — Mb for comparision. Errors are from the chiral fit only 
and the lattice data is adjusted for missing electromagnetism. 



for the Bd and B s . So to compare with experiment we 
shift A Bphy5 by -1 MeV to give 

M Bs - M Bl = 84(2) MeV 

in good agreement with experiment of 87.4(3) MeV 
(within 2 sigma). Reconstructing M B using our value 
for M Bs in Sec. [mA| gives M B = 5.283(2) (8) GeV The 



first error is from the chiral fit and the second is the error 
on M B with the detailed breakdown as in Sec. IIII A 



E. The D meson 

Our analysis of the D meson follows the same method 
as the B in the previous section. The splitting Md b — Mo 
is taken from a combined fit to all three charmed mesons, 
the results are given in Table |VIII| Systematic errors 
should be small in the splitting since the only difference 
between the states is the light quark mass, however we 
still see some lattice spacing dependence coming from 
the charm quark discretisation errors. We use the SU(2) 
HMxPT formula Eq. 16 with discretisation terms, this 
time including higher powers of a and with a scale set by 
m c 

(1.0+di(m c a) 2 +d 2 (m c a) i +d 3 (m c a) 6 +d 4 (m c a) s ). (18) 

Priors for g and the leading term are the same as above 
but priors for discretisation errors are 0.00(6) for d\ and 
0.0(2) for other <ij terms as in Sec. IIIB Since we do 



not have vector meson masses the experimental value 
140 MeV [22] is used for the hyperfine term in the fit 
function. The HISQ action has previously been shown 
to give results for hyperfine splittings in agreement with 
experiment [7]- 
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The result of the fit is shown in Fig. 10 including an 
adjustment for electromagnetism. The shift in the D s is 
1.3(7) MeV and the shifts in the D and D ± are -0.4 MeV 
and +1.3 MeV, which results in a total shift of 0.9 MeV 
in M Dl = (M Do + M D ±)/2. As for the B, tightening the 
prior on g to 0.1 also gives a consistent result, but in this 
case discretisation errors are significant so removing the 
di terms leads to a poorer fit. 

Our final result for the splitting is 

M Ds -M Dl = 101(3) MeV 

in agreement with the experimental splitting of 101.3(3) 
MeV [22] • When combined with our result for M& 3 above 
this gives Mn l — 1.869(3)(3)GeV, the first error being 
the chiral fitting error and the second the full error from 
Md 3 - 

Fig. [To] shows Mb s — Mb and M d s — Mr> on the same 
plot. It is clear that lattice QCD can distinguish the 
difference between these two small splittings. In HQET 
language it arises from the difference in the kinetic energy 
of the heavy quark in a heavy-strange meson compared 
to that in a heavy-light meson. We would expect this 
difference to be positive and contribute a larger amount 
for c quarks than b quarks, consistent with the increase 
seen. It is clear that lattice QCD successfully reproduces 
this effect. 



IV. HYPERFINE SPLITTINGS 

The hyperfine splitting between the ground state vec- 
tor and pseudo-scalar states is a particularly good test of 
a spectrum calculation. For heavy-light mesons this split- 
ting is proportional to the term -^^o ■ B in the NRQCD 
action so gives a direct check of the radiative corrections 
to C4. This is in contrast to the case in heavyonium 
where the hyperfine splitting is proportional to cf. The 
splitting also depends on higher order operators but in 
heavy-light systems these terms will be very small, unlike 
in bottomonium where the v 6 terms could be 10%. The 
splitting A^ p = Mb* — Ms q is very precise and including 
both the vector and pseudo-scalars in the same fit takes 
account of the correlations between the two. 

The bottomonium hyperfine splittings were calculated 
using our improved action in [3J, where we obtained 70(9) 
MeV for the IS hyperfine splitting and 0.499(42) for the 
ratio of the 2S and IS hyperfine splittings (which agreed 
well with subsequent experiment |23jV The error in both 



cases was dominated by the missing v terms. The heavy- 
light hyperfine splittings have previously been studied in 
Ref. [30] by considering ratios that were independent of 
c 4 . This resulted in a prediction of 53(7) MeV for the 
B c hyperfine splitting. The advantage of our current cal- 
culation is that the coefficients have been obtained by 
matching NRQCD to QCD at one loop, allowing the hy- 
perfine splittings to be determined directly without losing 
predictive power. Using the same action and C4 for both 



the bottomonium and B-meson calculations also allow us 
to make very different, independent checks. 

The results for the hyperfine splittings are given in 
lattice units in Table |XTV] for the Bi,B s and B c . Be- 
fore fitting the data we make a small correction for the b 
quark mass mistuning on each ensemble. The splittings 
are very insensitive to the light quark mass so retuning for 
m Sl m c will be negligible compared to other errors. The 
retuning assumes that the hyperfine splitting is inversely 
proportional to the b quark mass and is applied multi- 
plicatively using the tuned b quark mass values m^ hys 
calculated in [3] and listed in Table XIV There are two 
sources of error in hys , coming from the lattice spacing 
and the determination of the bottomonium kinetic mass 
values. Since a change in the lattice spacing would result 
in a change in the quark mass, the lattice spacing un- 
certainty is correlated with the scale uncertainty in the 
hyperfine splitting itself. To account for this correlation, 
we apply twice the lattice spacing error to the hyperfine 
splitting rather than adding them separately. The retun- 
ing factors are all less than 2% and are given in Table 

The dominant source of uncertainty in A"^ p is still the 
higher order correction to C4 which is now O(oq). To 
allow for this we apply a correlated systematic error to all 
the data points of size a 2 s where we take a s at a scale ir/a. 
Values for a s are: 0.275 on very coarse, 0.255 on coarse 
and 0.225 on fine [3]. The B s and Bi hyperfine splittings 
are fit to the same form as (Eq. |5])), allowing for 
lattice spacing, sea quark mass and cutoff dependence. 
The B c hyperfine is fit to the form used for As c in the 



hs method (Eq. (10)) in which we include discretisation 
errors with a scale set by m c . Priors are the same in all 
cases except for the prior on A^ p which is 0.5(5). 

The data and fit results are plotted in Fig. [Tljwith the 
data points adjusted for b quark mass mistuning but not 
including the correlated systematic error from C4. There 
is no noticable sea quark mass dependence and very little 
a dependence except for the B c case where the discretisa- 
tion errors come from the charm quark. The dependence 
of A"^ p on the light valence quark is also very small and 
not statistically significant. The results for the physical 
values from the fits are: 

A£ p = 50(3)MeV 



A£ s p = 52(3)MeV 
A h ™ = 60(6)MeV. 



(19) 

The other remaining source of error is the effect of v 6 
terms which are very small here. The full error budget 



for each splitting is given in Table XV Comparison to 
experiment (45.8(4) MeV [22]) for Bi shows good agree- 
ment. For B s the experimental results are not as ac- 
curate. In Fig. [TT] we use the experimental average of 
46.1(1.5) MeV [22 . This agrees with our value within 
2cr. 

The dominant error in the hyperfine splittings is still 
the uncertainty in the C4 coefficient which is reduced in 
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Set 
T~ 
2 
3 
4 
5 



0.0375(12) 

0.0377(9) 

0.0324(4) 

0.0309(4) 

0.0212(11) 



aA s s P 

0.03892(40) 

0.03705(47) 

0.03177(18) 

0.03102(16) 

0.02310(14) 



aA s c P 
0.03659(17) 

0.03658(13) 

0.03140(3) 

0.03137(2) 

0.02470(2) 



phvs 



3.297(11) (35) (7) (16) 
3.263(7)(35)(4)(16) 
2.696(4) (22) (7) (13) 
2.623(7) (22) (7) (13) 
1.893(6)(12)(5)(9) 



tuning 

1.000(5) 

1.000(5) 

0.987(5) 

0.999(6) 

1.009(5) 



TABLE XIV. Results for the hyperfine splittings Ag™ in lattice units for each ensemble, errors are statistical only. Column 5 
gives the tuned b quark masses calculated in [3] where the first two errors are from statistical and systematic errors respectively 
in the lattice spacing determination. The third and fourth errors are the statistical and systematic errors in the determining the 
Upsilon kinetic mass used for tuning am,b. The final column gives the multiplicative factor applied to each hyperfine splitting 
due to b quark mass mistuning. 
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FIG. 11. Results for the B, B s and B c hyperfine splittings. The data points include statistical, tuning and lattice spacing 
errors, taken as double the naive error as discussed in the text. The correlated a 2 error is not included on the data points. The 
result of the fit in each case is shown as a grey band and, where available, experimental values are given as black solid circles. 



TABLE XV. The full error budget for the hyperfine splittings, 
giving each error as a percentage of the final answer. The fit 
value is obtained including the statistical, scale and a 2 errors 
and their separate contribution to the error budget is distin- 
guished by fitting with and without the a 2 error, v 6 errors 
are included multiplicatively using the estimates in Sec. [TTJ 
The error from a, m q ^ sea , and arrib dependence is estimated 
from the fit. The error from mb tuning is estimated by fitting 
with and without the error on the tuning in Table |X"IV| 
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this calculation to 0(a 2 s ). Taking ratios of hyperfine 
splittings eliminates this error and also cancels errors 
from the lattice spacing and mistuning of the b quark 
mass. The remaining errors will be from missing v 6 terms 



of the B and B c hyperfine splittings to that of the B s . 
The fit function for A^P / A^ yp is the same as for the B^ 

hyperfine splitting and the fit function for A g yp /A } B yp is 
the same as for the B c hyperfine splitting. Priors on the 
ratios are taken to be 1.0(5). 
The results of the fits are 



Rb 



A]P_ 



0.993(33)(5) 



A 



hyp 



R Bo = -7& = 1.166(56)(23). 



(20) 



The first error is from statistics/fitting and the second 
is the systematic error that is dominated by missing v 6 
terms in the action. We take half the estimated size of 
v 6 terms as there should be some cancellation between 
the splittings. The full error budget is given in Table 



which are very small. Fig. 12 shows results for the ratio 



XV Our results are now precise enough that we are able 
to resolve the difference from 1.0 in the charm/strange 
hyperfine ratio, at the same time confirming our previous 
result [30 that this ratio is not far from 1. Our value 
for the ratio of light to strange hyperfine splittings is 1 
with an accuracy of 3% (equivalent to 1.5 MeV for this 
splitting). The experimental ratio, using the B s average 
above is 0.993(34). 
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The accuracy of the ratios above means that we can 
give an improved prediction of Mb* — Mg c . Multiplying 
the experimental average for B s by the ratio above gives: 



A^ p c = 54(3)MeV. 



(21) 



We take this as our final predicted value. Note that this 
is smaller than either the bottomonium or charmonium 
ground-state hyperfine splittings. 



V. AXIAL VECTOR AND SCALAR B c MESONS 

As discussed in Sec. [TTJ our £?-meson correlators con- 
tain oscillating terms corresponding to states of opposite 
parity. Hence our pseudoscalar correlators contain both 
0~ and + states, and the vector correlators contain l - 
and 1 + states. By using the fit form given in Eq. Q we 
can then extract the energies of these scalar and axial 
vector states from our fits. 

The splittings 



«*°b~°~ = aEB h - aE Bq 



(22) 
(23) 



are given in Tables |Vj [X] and XIII for the three mesons 
with q = s, c, I respectively. We find that, for the B, both 
states are above threshold for decay into Bit and for the 
B s the states are very close to threshold for BK decay, 
as was found in [9]. Since we do not have enough data 
to accurately estimate threshold effects in these cases we 
do not analyse them further, but they are included for 
completeness. 



The B c states, however, are far enough below thresh- 
old for decay to BD that we can reliably predict their 
masses. The remaining problem comes from identify- 
ing which states our results correspond to. From heavy 
quark spin symmetry, the "P-wave" heavy-light mesons 
come in two doublets, a + , 1 + pair coming from a light 
quark spin of ji — 1/2 and a 1 + , 2 + pair from ji — 3/2. 
Identifying our scalar state with the physical + state 
is unambiguous but the situation is not as clear for the 
axial-vector. Naively one would expect that we have cal- 
culated the lighter of the two states but without includ- 
ing a larger basis of operators this cannot be shown for 
certain. 

The results on coarse and fine ensembles for A° B ~° 



and A, 



are shown in Figs. 



13 



and [l4] respectively. 



form as in Eq. ( 11 ) 



The results are fit using the same 
with a prior of 0.5(5) on the physical value and the same 
priors as before for other parameters. As in the case of 
the radially excited states, we estimate that errors from 
missing relativistic corrections to the NRQCD action will 
be 1 MeV and that other systematic errors will be negli- 
gible. Our results for the splittings are then: 



B,. 

+ - 



= 429(13)(1) MeV 



(24) 



^j, = 410(13) (1) MeV, (25) 
where the first error is from the fit and the second is from 



NRQCD systematics. 



VI. DISCUSSION 



The results obtained here agree well with existing ex- 
periment and set improved levels of accuracy from a lat- 
tice QCD calculation. 



It is important to compare to other lattice QCD calcu- 
lations as well as to experiment because different lattice 
QCD methods have different systematic errors, partic- 
ularly if they use a different formalism for the quarks. 
Agreement then gives improved confidence in the error 
analysis. In Figs. 15 and 16 we compare existing re- 

l B r mesons from lattice 



suits for the masses of the B s am 
QCD, in which the quark masses are fixed from bottomo- 
nium, the r] c and the rj s . The comparison includes results 
from two very different formalisms for the b quark: the 
NRQCD formalism used here and in [9] and the HISQ 
formalism in which an extrapolation up to the b quark 
mass is made from lighter masses on lattices with a range 
of lattice spacings [lOl [IT] . The agreement between the 
different methods is good, within their total errors of 
around 10 MeV. 



VII. CONCLUSIONS 

We have presented results for the B meson spectrum 
using a perturbatively improved NRQCD action, very 
high statistics and gluon field configurations with an im- 
proved gluon action and including 2+1+1 flavours of 
HISQ sea quarks. We have improved upon and extended 
the previous results in Ref. [9] and, combined with our 
study of the Upsilon spectrum in Ref [2j , we have shown 
that our improved action gives accurate meson masses 
across a wide range of heavy mesons. Where we can 
compare, we see no significant differences with the re- 
sults of (5], so that the inclusion of c quarks in the sea 
has not produced any noticeable changes. 

The strongest improvement from our reduced system- 
atic errors can be seen in the hyperfine splittings which 
were previously dominated by missing radiative correc- 
tions. Our errors are now 3-6 MeV, giving an even more 
stringent test against experiment than for the bottomo- 
nium hyperfine splitting. The high statistics used in our 
calculation (32k correlators with 3 quark smearings) al- 
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FIG. 12. Ratios of B-meson hyperfine splittings. The result of the fit in each case is shown as a grey band and, where available, 
experimental values are given as black solid circles. 
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FIG. 13. Fit and lattice data on the fine and coarse ensembles 
for the splitting A 1 ^ _0 . Errors include statistics and scale 
uncertainty only. 



FIG. 14. Fit and lattice data on the fine and coarse ensembles 
for the splitting _1 . Errors include statistics and scale 
uncertainty only. 



lowed for the lightest B states to be reliably extracted 
and, with the light sea quark masses now available, con- 
sistent results were obtained for a range of reasonable 
chiral fit functions. This demonstration is particularly 
important for future determinations of f B which are cur- 
rently underway including ensembles with physical light 
quark masses. The calculation showed that lattice QCD 
could successfully resolve the change in splitting between 
heavy-strange and heavy-light meson masses as the quark 
mass is increased from c to b. The statistical precision 
of our correlators also allowed us to make the first QCD 
prediction of the radially excited B c states and two of 
the "P-wave" states. 

An overview of our results for the B-meson spectrum 
is shown in Fig. [T7] including the full error on each point. 
We find excellent agreement with the experimentally 



known pseudoscalar and vector states. In summary our 
results are: M Ba - M B , = 84(2) MeV, M Bs = 5.366(8) 
GeV, M Bc = 6.278(9) GeV, M Ds = 1.9697(33) GeV, and 
M Da - M D = 101(3) MeV. Our results for the B meson 
hyperfine splittings are M B » — M B — 50(3) MeV and 
M B , - M Bs = 52(3) MeV and we predict M B * - M Bc = 
54(3) MeV. Combining our results for the B c and the 
pseudoscalar radial splitting, we predict the mass of the 
B' c to be M B i = 6.894(19) Btat (8) syat GeV. Combining the 
B c , the hyperfine splitting and the vector radial splitting, 
we predict M B ,> = 6.922(19) atat (8) 3y5t GeV. Our predic- 
tion for the O+'state is M B « o = 6.707(14) stat (8) 3y3t GeV. 

Finally, in Fig. [18] we update the complete spectrum 
plot for gold-plated mesons to include the new results 
from this paper, as well as updated experimental val- 
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FIG. 15. A comparison of results for the B s meson mass from 
different formalisms for the 6 quark in lattice QCD. The exper- 
imental average value is given at the top with accompanying 
vertical lines. 
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FIG. 16. A comparison of results for the B c meson mass from 
different formalisms for the b quark in lattice QCD. In each 
case the result from the hh method is given above the result 
for the hs method. The experimental average value is given at 
the top with accompanying vertical lines. 
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FIG. 17. The spectrum of gold-plated B, B s and B c meson 
states from this calculation, compared to experiment where 
results exist. Predictions are marked with open red circles. 
None of the meson masses included here were used to tune 
parameters of the action so all the masses are parameter-free 
results from lattice QCD. 



ues. This plot summarises the coverage and the predic- 
tive power of lattice QCD calculations. 
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